home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 1.iso / toolbox / src / tutorials / OGLT / Examples / TerrainFollowing / src.hansen / culling.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-11-11  |  40.8 KB  |  1,418 lines
  1.  /**************************************************************************
  2.   *                                       *
  3.   *   Copyright (C) 1996, Silicon Graphics, Inc.               *
  4.   *              All Rights Reserved.                   *
  5.   *                                       *
  6.   *  The files in this subtree  contain  UNPUBLISHED  PROPRIETARY  SOURCE  *
  7.   *  CODE of Silicon Graphics, Inc.;  the  contents  of  these  files may  *
  8.   *  not be disclosed to third  parties,  copied  or  duplicated  in  any  *
  9.   *  form, in whole or in part,  without  the  prior  written  permission  * 
  10.   *  of  Silicon Graphics, Inc.                           *
  11.   *                                       *
  12.   *  RESTRICTED RIGHTS LEGEND:                           *
  13.   *  Use,  duplication  or  disclosure  by  the  Government is subject to  *
  14.   *  restrictions as set forth in subdivision (c)(1)(ii) of the Rights in  *
  15.   *  Technical Data and Computer Software clause at  DFARS  252.227-7013,  *
  16.   *  and/or in similar or successor  clauses in the FAR,  DOD or NASA FAR  *
  17.   *  Supplement.  Unpublished - rights reserved  under the Copyright Laws  *
  18.   *  of the United States.                           *
  19.   *                                       *
  20.   *  THIS SOFTWARE IS  PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND,  *
  21.   *  EXPRESS,  IMPLIED OR  OTHERWISE,  INCLUDING  WITHOUT LIMITATION, ANY  *
  22.   *  WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.      *
  23.   *                                       *
  24.   *  IN  NO  EVENT  SHALL  SILICON  GRAPHICS  BE  LIABLE FOR ANY SPECIAL,  *
  25.   *  INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF  ANY  KIND,  OR ANY  *
  26.   *  DAMAGES WHATSOEVER RESULTING FROM LOSS  OF  USE,  DATA  OR  PROFITS,  *
  27.   *  WHETHER OR NOT ADVISED OF  THE  POSSIBILITY  OF  DAMAGE,  AND ON ANY  *
  28.   *  THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR  *
  29.   *  PERFORMANCE OF THIS SOFTWARE.                       *
  30.   **************************************************************************/
  31.  
  32. /*
  33.     Author : Patrick Bouchaud    & Paul Hansen
  34.          galaad@neu.sgi.com    hansen@engr.sgi.com
  35. */
  36.  
  37. #include <stdlib.h>
  38. #include <stdio.h>
  39. #include <GL/gl.h>
  40. #include <bstring.h>
  41. #include <math.h>
  42. #include "data.h"
  43.  
  44. #include "culling.h"
  45.  
  46. #define ABS(x) ((x)<0?-(x):(x))
  47.  
  48. Segment* alloc_Segment(Segment*);
  49. Vertex* alloc_Vertex(Vertex*);
  50. Face* alloc_Face(Face*);
  51. Edge* alloc_Edge(Edge*);
  52. Volume* alloc_Volume(Volume*);
  53.  
  54. static GLvoid GetLOD(float, float, int*, int*);
  55.  
  56. /* These functions can be expanded to match image to terrain */
  57. #define LAT_CONV ((GLdouble)NUMTILE_S/NUM_LATITUDE/D_LATITUDE)
  58. #define LON_CONV ((GLdouble)NUMTILE_T/NUM_LONGITUDE/D_LONGITUDE)
  59. #define TILE_TO_ELEV_LAT(x) ((x)*(NUM_LATITUDE-1)*6144/6000/NUMTILE_S)
  60. #define TILE_TO_ELEV_LON(x) ((x)*(NUM_LONGITUDE-1)*6144/6000/NUMTILE_T)
  61.  
  62. extern float obsLat, obsLon, obsAlt, texture_select;
  63.  
  64. float cvec[2][2][3] = {{{1, 0, 0}, {0, 1, 0}}, {{0, 0, 1}, {1, 1, 0}}};
  65. /* float cvec[9][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 1, 0},
  66.             {0, 1, 1}, {1, 0, 1}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; */
  67. float row_bottom[NUMTILE_S][2];
  68. float row_top[NUMTILE_S][2];
  69. float row_perc_b[NUMTILE_S][2];
  70. float row_perc_t[NUMTILE_S][2];
  71. float col_left[NUMTILE_T][2];
  72. float col_right[NUMTILE_T][2];
  73. float col_perc_l[NUMTILE_T][2];
  74. float col_perc_r[NUMTILE_T][2];
  75. GLdouble combMat[16];
  76.  
  77. /*---------------------------------------------------------------------------* 
  78.  *    InitializeRowData -- establishes tile boundaries in lat/lon space      *
  79.  *---------------------------------------------------------------------------*/
  80.  
  81. void
  82. InitializeRowData(lat_stride, lon_stride)
  83. {
  84.     register int j, tmp;
  85.     float inv_latstride = (float)1.0 / lat_stride;
  86.     float inv_lonstride = (float)1.0 / lon_stride;
  87.  
  88.     for(j=0; j<NUMTILE_S; ++j) {
  89.     row_bottom[j][0] = TILE_TO_ELEV_LAT((float)j);
  90.     row_top[j][0] = row_bottom[j][1] = TILE_TO_ELEV_LAT((float)j+0.5);
  91.     row_top[j][1] = TILE_TO_ELEV_LAT((float)(j+1));
  92.  
  93.     tmp = (int)(row_bottom[j][0] * inv_latstride) * lat_stride;
  94.     if((float)tmp != row_bottom[j][0]) tmp += lat_stride;
  95.     row_perc_b[j][0] = ((float)tmp - row_bottom[j][0]) * inv_latstride;
  96.  
  97.     tmp = (int)(row_top[j][0] * inv_latstride) * lat_stride;
  98.     row_perc_t[j][0] = (row_top[j][0] - (float)tmp) * inv_latstride;
  99.  
  100.     tmp = (int)(row_bottom[j][1] * inv_latstride) * lat_stride;
  101.     if((float)tmp != row_bottom[j][1]) tmp += lat_stride;
  102.     row_perc_b[j][1] = ((float)tmp - row_bottom[j][1]) * inv_latstride;
  103.  
  104.     tmp = (int)(row_top[j][1] * inv_latstride) * lat_stride;
  105.     row_perc_t[j][1] = (row_top[j][1] - (float)tmp) * inv_latstride;
  106.     }
  107.     for(j=0; j<NUMTILE_T; ++j) {
  108.     col_left[j][0] = TILE_TO_ELEV_LON((float)j);
  109.     col_right[j][0] = col_left[j][1] = TILE_TO_ELEV_LON((float)j+0.5);
  110.     col_right[j][1] = TILE_TO_ELEV_LON((float)(j+1));
  111.  
  112.     tmp = (int)(col_left[j][0] * inv_lonstride) * lon_stride;
  113.     if((float)tmp != col_left[j][0]) tmp += lon_stride;
  114.     col_perc_l[j][0] = ((float)tmp - col_left[j][0]) * inv_lonstride;
  115.  
  116.     tmp = (int)(col_right[j][0] * inv_lonstride) * lon_stride;
  117.     col_perc_r[j][0] = (col_right[j][0] - (float)tmp) * inv_lonstride;
  118.  
  119.     tmp = (int)(col_left[j][1] * inv_lonstride) * lon_stride;
  120.     if((float)tmp != col_left[j][1]) tmp += lon_stride;
  121.     col_perc_l[j][1] = ((float)tmp - col_left[j][1]) * inv_lonstride;
  122.  
  123.     tmp = (int)(col_right[j][1] * inv_lonstride) * lon_stride;
  124.     col_perc_r[j][1] = (col_right[j][1] - (float)tmp) * inv_lonstride;
  125.     }
  126. }
  127.  
  128. /*---------------------------------------------------------------------------*
  129.  *    XCullDraw -- terrain is rendered as a series of textured tiles. Each tile
  130.  *     is rendered as 9 areas, some of which may not be in the viewing area.
  131.  *     The reader is referred to the "tiling_doc.sc" for diagrams and additional
  132.  *     description.
  133.  *---------------------------------------------------------------------------*/
  134.  
  135. #define LAT_METERS_TO_TILES(x) ((x)*NUMTILE_S*6000/6144/(NUM_LATITUDE-1)/D_LATITUDE)
  136. #define LON_METERS_TO_TILES(x) ((x)*NUMTILE_T*6000/6144/(NUM_LONGITUDE-1)/D_LONGITUDE)
  137. #define SCONV(x) (((x)-s_offset)*ISCALE_S)
  138. #define TCONV(x) (((x)-t_offset)*ISCALE_T)
  139.  
  140. #ifdef TEXTURED
  141. #define DO_VERTEX(lon, alt, lat) \
  142.     glTexCoord2f(SCONV(lat), TCONV(lon)); \
  143.     glVertex3f(lon*D_LONGITUDE, alt, lat*D_LATITUDE);
  144. #else
  145. #define DO_VERTEX(lon, alt, lat) \
  146.     glVertex3f(lon*D_LONGITUDE, alt, lat*D_LATITUDE);
  147. #endif
  148.  
  149. #define STRIP_SETUP                                \
  150.     bypass=0;                                    \
  151.     lon1 = l_lon[qrow][idx-1];                            \
  152.     if(lon1 >= right && lon2 >= right) bypass=1;                \
  153.     if(lon1 >= r_lon[qrow][idx-1] && lon2 >= r_lon[qrow][idx]) bypass=1;    \
  154.     if(!bypass) {                                \
  155.     lon = lon1 < lon2 ? lon1 : lon2;                    \
  156.     l_lon_i = (int)(lon * inv_lonstride) * lon_stride;            \
  157.     if(lon == col_left[tilecol][qcol])                    \
  158.         perc_l = col_perc_l[tilecol][qcol];                    \
  159.     else perc_l = 0.0;                            \
  160.     lon2 = r_lon[qrow][idx];                        \
  161.     lon1 = r_lon[qrow][idx-1];                        \
  162.     lon = lon1 > lon2 ? lon1 : lon2;                    \
  163.     if(lon > right)    {                            \
  164.         lon = right;                            \
  165.         perc_r = col_perc_r[tilecol][qcol];                    \
  166.         }                                    \
  167.     else perc_r = 0.0;                            \
  168.     r_lon_i = (int)(lon * inv_lonstride) * lon_stride;            \
  169.     if((float)r_lon_i != lon)                        \
  170.         r_lon_i += lon_stride;                        \
  171.                                         \
  172.     lon1 = (float)l_lon_i;                            \
  173.     lon2 = lon1 + lon_stride;                        \
  174.     alt3 = &Terrain[(int)lat1][(int)lon2];                    \
  175.     alt4 = &Terrain[(int)lat2][(int)lon2];                    \
  176.     alt1 = alt3 - lon_stride;                        \
  177.     alt2 = alt4 - lon_stride;                        \
  178.     }
  179.  
  180. int
  181. XCullDraw(
  182.     float altmin, 
  183.     float altmax, 
  184.     int lat_stride, 
  185.     int lon_stride
  186. )
  187. {
  188.     register float lat, lat1, lat2, lat3, lon, lon1, lon2, lon3, s_offset, t_offset;
  189.     int bypass, idx, j, tilerow, tilecol, qrow, qcol, minlod, maxlod;
  190.     float dlat_stride, lonmin, lonmax, latmin, latmax, unused;
  191.     float dlon_stride, inv_latstride, inv_lonstride;
  192.     float jminf, jmin_tile, jmaxf, jmax_tile, jval, iminf, imaxf;
  193.     static float l_lon[2][MAXGRID_LAT], r_lon[2][MAXGRID_LAT];
  194.     float lon_min[2], lon_max[2], right, dist, size;
  195.     int b_lat_i[2], t_lat_i, l_lon_i, r_lon_i;
  196.     float perc_l, perc_r, perc_t[2], perc_b[2], *alt1, *alt2, *alt3, *alt4, atmp, asave;
  197.     GLdouble vec1[4];
  198.  
  199.     xdim = (float)((NUM_LONGITUDE-1)/lon_stride*lon_stride)*D_LONGITUDE;
  200.     zdim = (float)((NUM_LATITUDE-1)/lat_stride*lat_stride)*D_LATITUDE;
  201.  
  202.     inv_latstride = (float)1.0 / lat_stride;
  203.     inv_lonstride = (float)1.0 / lon_stride;
  204.  
  205. /*
  206.     compute the 2 concave polygon resulting from the projection
  207.     of the viewing frustum onto the planes defined by Y=altmax, Y=altmin
  208. */
  209.     deleteVolume( pyramid );
  210.     pyramid = newFrustumPyramid();
  211.     BBoxClip( pyramid, 0., xdim, altmin, altmax, 0., zdim );
  212.  
  213. /*
  214.     get min and max x-values of the projected viewing frustrums
  215.     derive min and max j-indices into the elevation array.
  216.     After this, for each tile row, we gather data on the
  217.     min and max values at the culling boundary for each
  218.     "strip" to be rendered.
  219. */
  220.     getMinMaxLat(pyramid, &latmin, &latmax);
  221.     jmin_tile = LAT_METERS_TO_TILES(latmin);
  222.     jmax_tile = LAT_METERS_TO_TILES(latmax);
  223.     jminf = latmin/D_LATITUDE;
  224.     jmaxf = latmax/D_LATITUDE;
  225.     *(unsigned long*)&unused = 0xbfbfbfbf;
  226.     dlat_stride = lat_stride * D_LATITUDE;
  227.     for(jval=jminf, tilerow=(int)jmin_tile; tilerow<=(int)jmax_tile; ++tilerow) {
  228.     lon_max[0] = lon_max[1] = (float)-1;
  229.     lon_min[0] = lon_min[1] = (float)NUM_LONGITUDE;
  230.     for(qrow=0; qrow<2; ++qrow) {
  231.         l_lon[qrow][0] = r_lon[qrow][0] = unused;
  232.         if(jval >= row_top[tilerow][qrow]) continue;
  233.         if(jval >= jmaxf) continue;
  234.         b_lat_i[qrow] = (int)(jval * inv_latstride) * lat_stride;
  235.         if(jval == row_bottom[tilerow][qrow])
  236.         perc_b[qrow] = row_perc_b[tilerow][qrow];
  237.         else perc_b[qrow] = 0.0;
  238.  
  239.         jval = row_top[tilerow][qrow];
  240.         if(jval > jmaxf) {
  241.         jval = jmaxf;
  242.         perc_t[qrow] = 0.0;
  243.         }
  244.         else perc_t[qrow] = row_perc_t[tilerow][qrow];
  245.         t_lat_i = (int)(jval * inv_latstride) * lat_stride;
  246.         if((float)t_lat_i != jval)
  247.         t_lat_i += lat_stride;
  248.  
  249.         lat = b_lat_i[qrow] * D_LATITUDE;
  250.         for(idx=0, j=b_lat_i[qrow]; j<=t_lat_i; j+=lat_stride, idx++, lat+=dlat_stride) {
  251.         if(getMinMaxLon(pyramid, lat, &l_lon[qrow][idx], &r_lon[qrow][idx])) {
  252.             l_lon[qrow][idx] *= INV_DLON;
  253.             r_lon[qrow][idx] *= INV_DLON;
  254.             if(l_lon[qrow][idx] < lon_min[qrow]) lon_min[qrow] = l_lon[qrow][idx];
  255.             if(r_lon[qrow][idx] > lon_max[qrow]) lon_max[qrow] = r_lon[qrow][idx];
  256.             }
  257.         else l_lon[qrow][idx] = r_lon[qrow][idx] = -1.0;
  258.         }
  259.         l_lon[qrow][idx] = r_lon[qrow][idx] = unused;
  260.         if(l_lon[qrow][0] == -1.0) {
  261.         l_lon[qrow][0] = l_lon[qrow][1];
  262.         r_lon[qrow][0] = r_lon[qrow][1];
  263.         }
  264.         if(idx > 1 && l_lon[qrow][idx-1] == -1) {
  265.         l_lon[qrow][idx-1] = l_lon[qrow][idx-2];
  266.         r_lon[qrow][idx-1] = r_lon[qrow][idx-2];
  267.         }
  268.         }
  269.  
  270.     lonmin = INF(lon_min[0], lon_min[1]);
  271.     lonmax = SUP(lon_max[0], lon_max[1]);
  272.     iminf = lonmin * NUMTILE_T * 6000 / 6144 / (NUM_LONGITUDE-1);
  273.     imaxf = lonmax * NUMTILE_T * 6000 / 6144 / (NUM_LONGITUDE-1);
  274.  
  275.     /* data collection finished; now do rendering for this tile row */
  276.     glColor3f(1.0, 1.0, 1.0);
  277.     for(tilecol=(int)iminf; tilecol <= (int)imaxf; ++tilecol) {
  278. #ifdef TEXTURED
  279.         GetLOD(row_top[tilerow][0], col_right[tilecol][0], &minlod, &maxlod);
  280.         if(texture_select) {
  281.         glBindTextureEXT(GL_TEXTURE_2D, tilecol*NUMTILE_S+tilerow+1);
  282.         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_BASE_LEVEL_SGIS, minlod);
  283.         glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_LOD_SGIS, (float)minlod);
  284.         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAX_LEVEL_SGIS, maxlod);
  285.         }
  286. #endif
  287.         for(qcol=0; qcol<2; ++qcol) {
  288.         for(qrow=0; qrow<2; ++qrow) {
  289.             if(l_lon[qrow][0] == unused) continue;
  290.             if(lon_min[qrow] >= (right=col_right[tilecol][qcol]) ||
  291.                     lon_min[qrow] >= lon_max[qrow]) {
  292.             lon_min[qrow] = right;
  293.             continue;
  294.             }
  295.  
  296. #ifdef TEXTURED
  297.             s_offset = row_bottom[tilerow][qrow];
  298.             t_offset = col_left[tilecol][qcol];
  299.             if(texture_select)
  300.             glTexParameteri(GL_TEXTURE_2D, GL_QUAD_TEXTURE_SELECT_SGIS, (qcol<<1)|qrow);
  301.             else {
  302.             glBindTextureEXT(GL_TEXTURE_2D, ((tilecol*NUMTILE_S+tilerow<<2)|(qcol<<1)|qrow)+1);
  303.             glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_BASE_LEVEL_SGIS, minlod);
  304.             glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_LOD_SGIS, (float)minlod);
  305.             glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAX_LEVEL_SGIS, maxlod);
  306.             }
  307. #else
  308.             glColor3fv(cvec[qrow][qcol]);
  309. #endif
  310.             idx = 1;
  311.             lat1 = (float)b_lat_i[qrow];
  312.             lat2 = lat1 + lat_stride;
  313.             if(perc_b[qrow]) { /* do bottom tile stitching */
  314.             if((lon2 = l_lon[qrow][idx]) == unused) continue;    
  315.             STRIP_SETUP
  316.             if(!bypass) {
  317.                 lat3 = row_bottom[tilerow][qrow];
  318.  
  319.                 if(perc_l) { /* AREA #1 -- see "tiling_doc.sc" */
  320.                 lon3 = col_left[tilecol][qcol];
  321.                 glBegin(GL_TRIANGLE_FAN);
  322.                 if(perc_l + perc_b[qrow] > 1.0) {
  323.                     lat = lat1+perc_l*lat_stride;
  324.                     DO_VERTEX(lon3, *alt3+perc_l*(*alt2-*alt3), lat);
  325.                     }
  326.                 DO_VERTEX(lon3, *alt4+perc_l*(*alt2-*alt4), lat2);
  327.                 DO_VERTEX(lon2, *alt4, lat2);
  328.                 DO_VERTEX(lon2, atmp= *alt4+perc_b[qrow]*(*alt3-*alt4), lat3);
  329.                 if(perc_l + perc_b[qrow] > 1.0) {
  330.                     lon = lon1+perc_b[qrow]*lon_stride;
  331.                     DO_VERTEX(lon, *alt2+perc_b[qrow]*(*alt3-*alt2), lat3);
  332.                     }
  333.                 atmp = atmp + perc_l * ((*alt2+perc_b[qrow]*(*alt1-*alt2)) - atmp);
  334.                 DO_VERTEX(lon3, atmp, lat3);
  335.                 glEnd();
  336.                 lon1 = lon2; lon2 += lon_stride;
  337.                 alt1 = alt3; alt3 += lon_stride;
  338.                 alt2 = alt4; alt4 += lon_stride;
  339.                 }
  340.                 lon3 = lon1 + perc_b[qrow] * lon_stride;
  341.                 asave = *alt2 + perc_b[qrow] * (*alt1-*alt2);
  342.  
  343.                 /* AREA #2 -- see "tiling_doc.sc" */
  344.                 for(; lon2<=(float)r_lon_i-perc_r; lon1=lon2, lon2+=lon_stride, lon3+=lon_stride, 
  345.                         alt1=alt3, asave=atmp, alt3+=lon_stride, alt2=alt4, alt4+=lon_stride) {
  346.                 glBegin(GL_TRIANGLE_FAN);
  347.                 DO_VERTEX(lon1, *alt2, lat2);
  348.                 DO_VERTEX(lon2, *alt4, lat2);
  349.                 DO_VERTEX(lon2, atmp= *alt4+perc_b[qrow]*(*alt3-*alt4), lat3);
  350.                 DO_VERTEX(lon3, *alt2+perc_b[qrow]*(*alt3-*alt2), lat3);
  351.                 DO_VERTEX(lon1, asave, lat3);
  352.                 glEnd();
  353.                 }
  354.  
  355.                 if(perc_r) { /* AREA #3 -- see "tiling_doc.sc" */
  356.                 lon3 = right;
  357.                 glBegin(GL_TRIANGLE_FAN);
  358.                 DO_VERTEX(lon1, *alt2, lat2);
  359.                 DO_VERTEX(lon3, atmp= *alt2+perc_r*(*alt4-*alt2), lat2);
  360.                 atmp = atmp + perc_b[qrow] * ((*alt1+perc_r*(*alt3-*alt1)) - atmp);
  361.                 if(perc_r > perc_b[qrow]) {
  362.                     DO_VERTEX(lon3, atmp, lat3);
  363.                     lon = lon1+perc_b[qrow]*lon_stride;
  364.                     DO_VERTEX(lon, *alt2+perc_b[qrow]*(*alt3-*alt2), lat3);
  365.                     }
  366.                 else {
  367.                     lat = lat2-perc_r*lat_stride;
  368.                     DO_VERTEX(lon3, *alt2+perc_r*(*alt3-*alt2), lat);
  369.                     DO_VERTEX(lon3, atmp, lat3);
  370.                     }
  371.                 DO_VERTEX(lon1, *alt2+perc_b[qrow]*(*alt1-*alt2), lat3);
  372.                 glEnd();
  373.                 }
  374.                 }
  375.             lat1 = lat2; lat2 += lat_stride;
  376.             if(right > l_lon[qrow][0]) 
  377.                 l_lon[qrow][0] = right;
  378.             idx++;
  379.             }
  380.             for(; ; idx++, lat1=lat2, lat2+=lat_stride) { /* now do main section */
  381.             if((lon2 = l_lon[qrow][idx]) == unused) break;    
  382.             if(l_lon[qrow][idx+1]==unused && perc_t[qrow]) break;
  383.             STRIP_SETUP
  384.             if(!bypass) {
  385.  
  386.                 if(perc_l) { /* AREA #4 -- see "tiling_doc.sc" */
  387.                 lon3 = col_left[tilecol][qcol];
  388.                 glBegin(GL_TRIANGLE_FAN);
  389.                 DO_VERTEX(lon2, *alt3, lat1);
  390.                 DO_VERTEX(lon3, *alt3+perc_l*(*alt1-*alt3), lat1);
  391.                 lat = lat1+perc_l*lat_stride;
  392.                 DO_VERTEX(lon3, *alt3+perc_l*(*alt2-*alt3), lat);
  393.                 DO_VERTEX(lon3, *alt4+perc_l*(*alt2-*alt4), lat2);
  394.                 DO_VERTEX(lon2, *alt4, lat2);
  395.                 glEnd();
  396.                 lon1 = lon2; lon2 += lon_stride;
  397.                 alt1 = alt3; alt3 += lon_stride;
  398.                 alt2 = alt4; alt4 += lon_stride;
  399.                 }
  400.  
  401.                 /* AREA #5 -- see "tiling_doc.sc" */
  402.                 glBegin(GL_TRIANGLE_STRIP);
  403.                 for(; lon1<=(float)r_lon_i-perc_r; lon1+=lon_stride, alt1+=lon_stride, alt2+=lon_stride) {
  404.                 DO_VERTEX(lon1, *alt1, lat1);
  405.                 DO_VERTEX(lon1, *alt2, lat2);
  406.                 }
  407.                 glEnd();
  408.  
  409.                 if(perc_r) { /* AREA #6 -- see "tiling_doc.sc" */ 
  410.                 lon2 = lon1; lon1 -= lon_stride;
  411.                 alt3 = alt1; alt1 -= lon_stride;
  412.                 alt4 = alt2; alt2 -= lon_stride;
  413.                 lon3 = right;
  414.                 glBegin(GL_TRIANGLE_FAN);
  415.                 DO_VERTEX(lon1, *alt2, lat2);
  416.                 DO_VERTEX(lon3, *alt2+perc_r*(*alt4-*alt2), lat2);
  417.                 lat = lat2+perc_r*(lat1-lat2);
  418.                 DO_VERTEX(lon3, *alt2+perc_r*(*alt3-*alt2), lat);
  419.                 DO_VERTEX(lon3, *alt1+perc_r*(*alt3-*alt1), lat1);
  420.                 DO_VERTEX(lon1, *alt1, lat1);
  421.                 glEnd();
  422.                 }
  423.                 }
  424.             if(right > l_lon[qrow][idx-1]) 
  425.                 l_lon[qrow][idx-1] = right;
  426.             }
  427.             if(perc_t[qrow]) { /* do top tile stitching */
  428.             if((lon2 = l_lon[qrow][idx]) == unused) continue;    
  429.             STRIP_SETUP
  430.             if(!bypass) {
  431.                 lat3 = row_top[tilerow][qrow];
  432.  
  433.                 if(perc_l) { /* AREA #7 -- see "tiling_doc.sc" */
  434.                 lon3 = col_left[tilecol][qcol];
  435.                 glBegin(GL_TRIANGLE_FAN);
  436.                 DO_VERTEX(lon2, *alt3, lat1);
  437.                 DO_VERTEX(lon3, atmp= *alt3+perc_l*(*alt1-*alt3), lat1);
  438.                 atmp = atmp + perc_t[qrow] * ((*alt4+perc_l*(*alt2-*alt4)) - atmp);
  439.                 if(perc_l > perc_t[qrow]) {
  440.                     DO_VERTEX(lon3, atmp, lat3);
  441.                     lon = lon2-perc_t[qrow]*lon_stride;
  442.                     DO_VERTEX(lon, *alt3+perc_t[qrow]*(*alt2-*alt3), lat3);
  443.                     }
  444.                 else {
  445.                     lat = lat1+perc_l*lat_stride;
  446.                     DO_VERTEX(lon3, *alt3+perc_l*(*alt2-*alt3), lat);
  447.                     DO_VERTEX(lon3, atmp, lat3);
  448.                     }
  449.                 DO_VERTEX(lon2, *alt3+perc_t[qrow]*(*alt4-*alt3), lat3);
  450.                 glEnd();
  451.                 lon1 = lon2; lon2 += lon_stride;
  452.                 alt1 = alt3; alt3 += lon_stride;
  453.                 alt2 = alt4; alt4 += lon_stride;
  454.                 }
  455.                 lon3 = lon2 - perc_t[qrow] * lon_stride;
  456.                 asave = *alt1 + perc_t[qrow] * (*alt2-*alt1);
  457.  
  458.                 /* AREA #8 -- see "tiling_doc.sc" */
  459.                 for(; lon2<=(float)r_lon_i-perc_r; lon1=lon2, lon2+=lon_stride, lon3+=lon_stride, 
  460.                             alt1=alt3, alt3+=lon_stride, alt2=alt4, alt4+=lon_stride) {
  461.                 glBegin(GL_TRIANGLE_FAN);
  462.                 DO_VERTEX(lon2, *alt3, lat1);
  463.                 DO_VERTEX(lon1, *alt1, lat1);
  464.                 DO_VERTEX(lon1, asave, lat3);
  465.                 DO_VERTEX(lon3, *alt3+perc_t[qrow]*(*alt2-*alt3), lat3);
  466.                 DO_VERTEX(lon2, asave= *alt3+perc_t[qrow]*(*alt4-*alt3), lat3);
  467.                 glEnd();
  468.                 }
  469.  
  470.                 if(perc_r) { /* AREA #9 -- see "tiling_doc.sc" */
  471.                 lon3 = col_right[tilecol][qcol];
  472.                 glBegin(GL_TRIANGLE_FAN);
  473.                 if(perc_r + perc_t[qrow] > 1.0) {
  474.                     lat = lat2-perc_r*lat_stride;
  475.                     DO_VERTEX(lon3, *alt2+perc_r*(*alt3-*alt2), lat);
  476.                     }
  477.                 DO_VERTEX(lon3, *alt1+perc_r*(*alt3-*alt1), lat1);
  478.                 DO_VERTEX(lon1, *alt1, lat1);
  479.                 DO_VERTEX(lon1, atmp= *alt1+perc_t[qrow]*(*alt2-*alt1), lat3);
  480.                 if(perc_r + perc_t[qrow] > 1.0) {
  481.                     lon = lon2-perc_t[qrow]*lon_stride;
  482.                     DO_VERTEX(lon, *alt3+perc_t[qrow]*(*alt2-*alt3), lat3);
  483.                     }
  484.                 atmp = atmp + perc_r * ((*alt3+perc_t[qrow]*(*alt4-*alt3)) - atmp);
  485.                 DO_VERTEX(lon3, atmp, lat3);
  486.                 glEnd();
  487.                 }
  488.                 }
  489.             if(right > l_lon[qrow][idx-1]) 
  490.                 l_lon[qrow][idx-1] = right;
  491.             idx++;
  492.             }
  493.             if(right > l_lon[qrow][idx])
  494.             l_lon[qrow][idx-1] = right;
  495.             }
  496.         }
  497.         }
  498.     /* fprintf(stderr, "one tile row\n"); */
  499.     }
  500. }
  501.  
  502. /*---------------------------------------------------------------------------* 
  503.  *    GetLOD -- calculates a range of mipmap levels to be used for tile 
  504.  *     Calculates distance to point (lat,lon), which is taken as a ratio
  505.  *     to a standard distance at which level 0 appears to be full size,
  506.  *     i.e., neither mag nor min. We then take the log of the ratio by
  507.  *     extracting the exponent of the floating point value.
  508.  *---------------------------------------------------------------------------*/
  509.  
  510. static GLvoid
  511. GetLOD(
  512.     float lat,
  513.     float lon,
  514.     int *minlod,
  515.     int *maxlod
  516.     ) {
  517.     register GLfloat atmp, dist;
  518.     GLfloat size;
  519.  
  520.     atmp = obsAlt/*  - Terrain[(int)lat][(int)lon] */;
  521.     if(atmp < 0) atmp = 0;
  522.     lat = ABS(lat*D_LATITUDE-obsLat);
  523.     lon = ABS(lon*D_LONGITUDE-obsLon);
  524.     if(lat > lon)
  525.     dist = sqrt(lat*lat+atmp*atmp);
  526.     else dist = sqrt(lon*lon+atmp*atmp);
  527.     if(dist==0) size = 256.0;
  528.     else size = dist/(7000 * atmp/dist);
  529.     *minlod = ((*(unsigned long*)&size>>23)&0xff)-127;
  530.     if(*minlod < 0) *minlod = 0;
  531.     *maxlod = *minlod+2;
  532.     if(*maxlod > 4) *maxlod = 8;
  533.     /*
  534.     fprintf(stderr, "lat: %f lon: %f atmp: %f size: %f minlod: %ld\n", lat, lon, obsAlt, size, *minlod);
  535.     */
  536. }
  537.  
  538. /*---------------------------------------------------------------------------*/
  539. /*        newFrustumPyramid
  540. /*---------------------------------------------------------------------------*/
  541.  
  542. #define    FrontLeftBottom    0
  543. #define    BackLeftBottom    1
  544. #define    FrontRightBottom    2
  545. #define    BackRightBottom    3
  546. #define    FrontRightTop    4
  547. #define    BackRightTop    5
  548. #define    FrontLeftTop    6
  549. #define    BackLeftTop        7
  550.  
  551. #define    LeftFace        0
  552. #define    RightFace        1
  553. #define    TopFace        2
  554. #define    BottomFace        3
  555. #define    FrontFace        4
  556. #define    BackFace        5
  557.  
  558. static Volume *
  559. newFrustumPyramid()
  560. {
  561.     GLint    i, vp[4];
  562.     GLdouble    projMat[16], modelMat[16],
  563.         x0, y0, z0, x1, y1, z1,
  564.         *normal;
  565.     Volume    *volume;
  566.     Face    *Polygon[6], *face;
  567.     Edge    *edge;
  568.     Segment    *segment;
  569.     Vertex    *Corner[8], *v, *v0, *v1, *v2, *v3;
  570.  
  571. /*
  572.     Allocate the data structure for the new volume :
  573.     ------------------------------------------------
  574.     pyramid = 8 vertex / 6 faces / 12 segments
  575. */
  576.     volume = alloc_Volume( NULL );
  577.  
  578.     volume->numVertex = 8;
  579.     for (i=0; i<8; i++) Corner[i] = alloc_Vertex(NULL);
  580.     v = volume->firstVertex = Corner[0];
  581.     v->last = NULL;
  582.     for (i=1; i<8; i++)
  583.     {
  584.     v->next = Corner[i];
  585.     v->next->last = v;
  586.     v = v->next;
  587.     }
  588.     v->next = NULL;
  589.     volume->lastVertex = v;
  590.  
  591.     volume->numSegment = 0;
  592.     volume->firstSegment = volume->lastSegment = NULL;
  593.  
  594.     volume->numFace = 6;
  595.     for (i=0; i<6; i++) Polygon[i] = newFace( 4, NULL );
  596.     face = volume->firstFace = Polygon[0];
  597.     face->last = NULL;
  598.     for (i=1; i<6; i++)
  599.     {
  600.     face->next = Polygon[i];
  601.     face->next->last = face;
  602.     face = face->next;
  603.     }
  604.     face->next = NULL;
  605.     volume->lastFace = face;
  606.  
  607. /*
  608.     compute the 3D coordinates of the 8 corners of the viewing
  609.     ----------------------------------------------------------
  610.     frustum, as resulting from the gluUnProject( current matrix )
  611.     of the corners of the current viewport, with zval={0.,1.}
  612. */
  613.     glGetDoublev( GL_PROJECTION_MATRIX, projMat );
  614.     glGetDoublev( GL_MODELVIEW_MATRIX, modelMat );
  615.     glGetIntegerv( GL_VIEWPORT, vp );
  616.  
  617.     x0 = (GLdouble) vp[0];
  618.     y0 = (GLdouble) vp[1];
  619.     x1 = (GLdouble) (vp[0]+vp[2]-1);
  620.     y1 = (GLdouble) (vp[1]+vp[3]-1);
  621.  
  622. #define setVertex( wx, wy, wz, angle ) \
  623.     v = Corner[angle];\
  624.     v->in = 1;\
  625.     gluUnProject( wx, wy, wz, modelMat, projMat, vp, &v->x, &v->y, &v->z );
  626.  
  627.     setVertex( x0, y0, 0., FrontLeftBottom );
  628.     setVertex( x0, y0, 1., BackLeftBottom );
  629.     setVertex( x1, y0, 0., FrontRightBottom );
  630.     setVertex( x1, y0, 1., BackRightBottom );
  631.     setVertex( x1, y1, 0., FrontRightTop );
  632.     setVertex( x1, y1, 1., BackRightTop );
  633.     setVertex( x0, y1, 0., FrontLeftTop );
  634.     setVertex( x0, y1, 1., BackLeftTop );
  635.  
  636. /*
  637.     now compute the 6 normals of each face of the viewing frustum
  638.     -------------------------------------------------------------
  639.     from the previously computed corners.
  640.     The normals are NOT normalized, and are computed according to 
  641.     normal = (corner0,corner1)^(corner0,corner3)
  642.  
  643.     the viewing frustum is defined by the 6 inequations :
  644.     normal[face].(corner[face],M) >= 0
  645.     we compute now the coefficients of these equations
  646. */
  647. #define setFace( side, angle0, angle1, angle2, angle3 )\
  648.     \
  649.     face = Polygon[side];\
  650.     v0 = Corner[angle0];\
  651.     v1 = Corner[angle1];\
  652.     v2 = Corner[angle2];\
  653.     v3 = Corner[angle3];\
  654.     x0 = v1->x-v0->x; y0 = v1->y-v0->y; z0 = v1->z-v0->z;\
  655.     x1 = v3->x-v0->x; y1 = v3->y-v0->y; z1 = v3->z-v0->z;\
  656.     normal = face->equation;\
  657.     normal[0] = y0*z1 - y1*z0;\
  658.     normal[1] = x1*z0 - x0*z1;\
  659.     normal[2] = x0*y1 - x1*y0;\
  660.     normal[3] = -(v0->x*normal[0]+v0->y*normal[1]+v0->z*normal[2]);\
  661.     \
  662.     edge = face->firstEdge;\
  663.     edge->segment = newSegment( volume, v0, v1 );\
  664.     edge->inverted = (edge->segment->start==v1);\
  665.     edge = edge->next;\
  666.     edge->segment = newSegment( volume, v1, v2 );\
  667.     edge->inverted = (edge->segment->start==v2);\
  668.     edge = edge->next;\
  669.     edge->segment = newSegment( volume, v2, v3 );\
  670.     edge->inverted = (edge->segment->start==v3);\
  671.     edge = edge->next;\
  672.     edge->segment = newSegment( volume, v3, v0 );\
  673.     edge->inverted = (edge->segment->start==v0);
  674.  
  675.     face = volume->firstFace;
  676.     setFace( LeftFace,
  677.     BackLeftBottom, BackLeftTop,
  678.     FrontLeftTop, FrontLeftBottom
  679.     );
  680.  
  681.     face = face->next;
  682.     setFace( RightFace,
  683.     BackRightBottom, FrontRightBottom,
  684.     FrontRightTop, BackRightTop
  685.     );
  686.  
  687.     face = face->next;
  688.     setFace( FrontFace,
  689.     FrontLeftBottom, FrontLeftTop,
  690.     FrontRightTop, FrontRightBottom
  691.     );
  692.  
  693.     face = face->next;
  694.     setFace( BottomFace,
  695.     BackLeftBottom, FrontLeftBottom,
  696.     FrontRightBottom, BackRightBottom
  697.     );
  698.  
  699.     face = face->next;
  700.     setFace( TopFace,
  701.     BackLeftTop, BackRightTop,
  702.     FrontRightTop, FrontLeftTop
  703.     );
  704.  
  705.     face = face->next;
  706.     setFace( BackFace,
  707.     BackLeftBottom, BackRightBottom,
  708.     BackRightTop, BackLeftTop
  709.     );
  710.  
  711.     return volume;
  712. }
  713.  
  714. /*---------------------------------------------------------------------------*/
  715. /*    BBoxClip
  716. /*---------------------------------------------------------------------------*/
  717. static void
  718. BBoxClip( Volume *volume,
  719.     float xmin, float xmax,
  720.     float ymin, float ymax,
  721.     float zmin, float zmax )
  722. {
  723.     GLdouble equation[4];
  724.  
  725. #if 1
  726.     equation[0] = 1.;    /* x-xmin >=0 */
  727.     equation[1] = 0.;
  728.     equation[2] = 0.;
  729.     equation[3] = -xmin;
  730.     ClipVolume( volume, equation );
  731.  
  732.     equation[0] = -1.;    /* -x+xmax >= 0 */
  733.     equation[3] = xmax;
  734.     ClipVolume( volume, equation );
  735.  
  736.     equation[0] = 0.;    /* y-ymin >= 0 */
  737.     equation[1] = 1.;
  738.     equation[3] = -ymin;
  739.     ClipVolume( volume, equation );
  740.  
  741.     equation[1] = -1.;    /* -y+ymax >= 0 */
  742.     equation[3] = ymax;
  743.     ClipVolume( volume, equation );
  744.  
  745.     equation[1] = 0.;    /* z-zmin >= 0 */
  746.     equation[2] = 1.;
  747.     equation[3] = -zmin;
  748.     ClipVolume( volume, equation );
  749.  
  750.     equation[2] = -1.;    /* -z+zmax >= 0 */
  751.     equation[3] = zmax;
  752.     ClipVolume( volume, equation );
  753. #else
  754.     equation[0] = 0.;    /* -y+ymax >= 0 */
  755.     equation[1] = -1.;
  756.     equation[2] = 0.;
  757.     equation[3] = ymax;
  758.     ClipVolume( volume, equation );
  759. #endif
  760. }
  761.  
  762. /*---------------------------------------------------------------------------*/
  763. /*    ClipLine3D
  764. /*---------------------------------------------------------------------------*/
  765. /*
  766.     compute the 3D intersection of the
  767.     line: (v1,v2)
  768.     line: eq[0]*v->x + eq[1]*v->y + eq[2]*v->z + eq[3] = 0
  769.     (purely mathematical)
  770. */
  771. static Vertex *
  772. ClipLine3D( Vertex *v1, Vertex *v2, GLdouble *eq )
  773. {
  774.     GLdouble t;
  775.     Vertex *v = alloc_Vertex( NULL );
  776.     v->in = 1;
  777.  
  778.     t = equate(v1,eq);
  779.     if ( t == 0. )
  780.     {
  781.     v->x = v1->x;
  782.     v->y = v1->y;
  783.     v->z = v1->z;
  784.     }
  785.     else
  786.     {
  787.     t = t/(t-equate(v2,eq));
  788.     v->x = v1->x + t*( v2->x - v1->x );
  789.     v->y = v1->y + t*( v2->y - v1->y );
  790.     v->z = v1->z + t*( v2->z - v1->z );
  791.     }
  792.     return v;
  793. }
  794.  
  795. /*---------------------------------------------------------------------------*/
  796. /*    ClipVolume
  797. /*---------------------------------------------------------------------------*/
  798. static void
  799. ClipVolume( Volume *volume, GLdouble *equation )
  800. {
  801.     int    i, j, numf, numv, nums;
  802.     Vertex    *v, *v0, *v1;
  803.     Segment    *segment;
  804.     Face    *face;
  805.     Edge    *edge, *edge0, *edge1;
  806.  
  807.     if (volume==NULL) return;
  808.     if ((!volume->numVertex) || (!volume->numSegment) || (!volume->numFace))
  809.     return;
  810.  
  811.     volume->numNewSegment = 0;
  812.     volume->firstNewSegment = NULL;
  813. /*
  814.     test every vertex against clipping plane
  815.     ----------------------------------------
  816. */
  817.     numv = volume->numVertex;
  818.     v = volume->firstVertex;
  819.     for (i=0; i<numv; i++, v=v->next)
  820.     v->in = (equate(v,equation) >= 0);
  821. /*
  822.     Now, rework segments
  823.     --------------------
  824. */
  825.     nums = volume->numSegment;
  826.     segment = volume->firstSegment;
  827.     for (i=0; i<nums; i++, segment=segment->next)
  828.     {
  829.     if (segment->start->in)
  830.         segment->state = (segment->end->in) ? CULLED_IN : CLIPPED_END;
  831.     else
  832.         segment->state = (segment->end->in) ? CLIPPED_START : CULLED_OUT;
  833.  
  834.     if ((segment->state==CLIPPED_START) || (segment->state==CLIPPED_END))
  835.     {
  836.         v = ClipLine3D( segment->start, segment->end, equation);
  837.         volume->lastVertex->next = v;
  838.         v->next = NULL;
  839.         v->last = volume->lastVertex;
  840.         volume->lastVertex = v;
  841.         volume->numVertex++;
  842.     
  843.         if (segment->state==CLIPPED_START) segment->start = v;
  844.         else segment->end = v;
  845.     }
  846.     }
  847.  
  848. /*
  849.     clip every face of the Volume
  850.     -----------------------------
  851. */
  852. #define    GETSTATE(edge)\
  853.     ((edge->inverted) ?\
  854.     ((edge->segment->state==CLIPPED_END) ?\
  855.         CLIPPED_START\
  856.         :\
  857.         ((edge->segment->state==CLIPPED_START) ?\
  858.         CLIPPED_END\
  859.         :\
  860.         edge->segment->state\
  861.         )\
  862.     )\
  863.     :\
  864.     edge->segment->state\
  865.     )
  866.  
  867.     numf = volume->numFace;
  868.     face = volume->firstFace;
  869.     for (i=0; i<numf; i++, face=face->next)
  870.     {
  871.     Edge    *edge = face->firstEdge;
  872.     int    nume = face->numEdge;
  873.     int    state = GETSTATE(edge);
  874.  
  875.     switch    (state)
  876.     {
  877.         case CULLED_IN:
  878.         case CLIPPED_END:
  879.         if (state!=CLIPPED_END) do
  880.         {
  881.             edge = edge->next;
  882.             state = GETSTATE(edge);
  883.         }
  884.         while ((state!=CLIPPED_END) && (edge!=face->firstEdge));
  885.         if ( state == CLIPPED_END )
  886.         {
  887.             edge0 = edge;
  888.             edge = edge->next;
  889.             state = GETSTATE(edge);
  890.             while ((state!=CLIPPED_START) && (edge!=edge0))
  891.             {
  892.             edge = edge->next;
  893.             state = GETSTATE(edge);
  894.             }
  895.             if (state != CLIPPED_START)
  896.             {
  897.             fprintf(stderr,
  898.                 "ERROR: ClipVolume: unable to clip face\n");
  899.             break;
  900.             }
  901.     
  902.             segment = edge0->segment;
  903.             v0 = (edge0->inverted) ? segment->start : segment->end;
  904.             segment = edge->segment;
  905.             v1 = (edge->inverted) ? segment->end : segment->start;
  906.             segment = newSegment( volume, v0, v1 );
  907.     
  908.             edge1 = alloc_Edge( NULL );
  909.             edge1->segment = segment;
  910.             edge1->inverted = (segment->start==v1);
  911.             edge1->next = edge;
  912.             edge1->last = edge0;
  913.             edge = edge->last;
  914.             while (edge != edge0)
  915.             {
  916.             edge = edge->last;
  917.             alloc_Edge( edge->next );
  918.             nume--;
  919.             }
  920.             edge1->next->last = edge1;
  921.             edge1->last->next = edge1;
  922.     
  923.             face->firstEdge = edge1;
  924.             face->numEdge = nume+1;
  925.         }
  926.         break;
  927.  
  928.         case CULLED_OUT:
  929.         case CLIPPED_START:
  930.         if (state != CLIPPED_START) do
  931.         {
  932.             edge = edge->next;
  933.             state = GETSTATE(edge);
  934.         }
  935.         while ((state!=CLIPPED_START) && (edge!=face->firstEdge));
  936.         if ( state != CLIPPED_START )
  937.         {
  938.             for (j=0; j<nume; j++)
  939.             {
  940.             Edge    *nextEdge = edge->next;
  941.             alloc_Edge( edge );
  942.             edge = nextEdge;
  943.             }
  944.             face->numEdge = 0;
  945.             face->firstEdge = NULL;
  946.         }
  947.         else
  948.         {
  949.             edge0 = edge;
  950.             edge = edge->next;
  951.             state = GETSTATE(edge);
  952.             while ((state != CLIPPED_END) && (edge!=edge0))
  953.             {
  954.             edge = edge->next;
  955.             state = GETSTATE(edge);
  956.             }
  957.             if (state != CLIPPED_END)
  958.             {
  959.             fprintf(stderr,
  960.                 "ERROR: ClipVolume: unable to clip face\n");
  961.             break;
  962.             }
  963.             segment = edge->segment;
  964.             v0 = (edge->inverted) ? segment->start : segment->end;
  965.             segment = edge0->segment;
  966.             v1 = (edge0->inverted) ? segment->end : segment->start;
  967.             segment = newSegment( volume, v0, v1 );
  968.  
  969.             edge1 = alloc_Edge( NULL );
  970.             edge1->segment = segment;
  971.             edge1->inverted = (segment->start==v1);
  972.             edge1->next = edge0;
  973.             edge1->last = edge;
  974.             edge = edge->next;
  975.             while (edge!=edge0)
  976.             {
  977.             edge = edge->next;
  978.             alloc_Edge( edge->last );
  979.             nume--;
  980.             }
  981.             edge1->next->last = edge1;
  982.             edge1->last->next = edge1;
  983.  
  984.             face->firstEdge = edge1;
  985.             face->numEdge = nume+1;
  986.         }
  987.         break;
  988.     }
  989.     }
  990. /*
  991.     build a new face, made of all the newly created segments
  992.     --------------------------------------------------------
  993. */
  994. #ifdef DEV
  995.     fprintf( stderr, "numNewSegment = %d\n", volume->numNewSegment );
  996. #endif
  997.     if ((nums=volume->numNewSegment) > 0)
  998.     {
  999.     face = newFace( nums, equation );
  1000.     edge = face->firstEdge;
  1001.     segment = volume->firstNewSegment;
  1002.     for (i=0; i<nums; i++)
  1003.     {
  1004.         edge->segment = segment;
  1005.         segment = segment->next;
  1006.         edge = edge->next;
  1007.     }
  1008. /*
  1009.     re-order the edges, so that each segment starts with
  1010.     the end of the previous segment (of the previous edge)
  1011. */
  1012.     edge0 = face->firstEdge;
  1013.     edge0->inverted = 0;
  1014.     for (i=0; i<nums-1; i++)
  1015.     {
  1016.         int found;
  1017.         segment = edge0->segment;
  1018.         v = (edge0->inverted) ? segment->start : segment->end;
  1019.         edge = edge0->next;
  1020.         for (found=0, j=i+1; (!found) && (j<nums); j++)
  1021.         {
  1022.         segment = edge->segment;
  1023.         if ((segment->start == v) || (segment->end == v))
  1024.         {
  1025.             found = 1;
  1026.             edge->inverted = (segment->end == v);
  1027.         }
  1028.         else edge = edge->next;
  1029.         }
  1030.         if (found)
  1031.         {
  1032.         edge->next->last = edge->last;
  1033.         edge->last->next = edge->next;
  1034.         edge->next = edge0->next;
  1035.         edge->last = edge0;
  1036.         edge->next->last = edge;
  1037.         edge->last->next = edge;
  1038.         edge0 = edge;
  1039.         }
  1040.         else
  1041.         {
  1042.         fprintf(stderr,
  1043.             "ERROR::ClipVolume:: unable to close new face\n");
  1044.         edge0 = edge0->next;
  1045.         }
  1046.     }
  1047.     face->last = volume->lastFace;
  1048.     face->next = NULL;
  1049.     volume->lastFace->next = face;
  1050.     volume->lastFace = face;
  1051.     volume->numFace++;
  1052.     }
  1053. /*
  1054.     Finally remove all the useless data :
  1055.     --------------------------------------
  1056.     vertex removed if !(vertex->in)
  1057.     segment removed if (segment->state==CULLED_OUT)
  1058.     face removed if (face->numEdge==0)
  1059. */
  1060. #define    CLEANUP( Type, Func, condition, firstItem, lastItem, numItem )\
  1061.     {\
  1062.     Type *ptr = firstItem;\
  1063.     int num = numItem;\
  1064.     for (i=0; i<num; i++)\
  1065.     {\
  1066.         Type *nextptr = ptr->next;\
  1067.         if (condition)\
  1068.         {\
  1069.         Type *lastptr = ptr->last;\
  1070.         if (lastptr) lastptr->next = nextptr;\
  1071.         if (nextptr) nextptr->last = lastptr;\
  1072.         if (firstItem==ptr) firstItem=nextptr;\
  1073.         if (lastItem==ptr) lastItem=lastptr;\
  1074.         Func( ptr );\
  1075.         numItem--;\
  1076.         }\
  1077.         ptr = nextptr;\
  1078.     }\
  1079.     }
  1080.  
  1081.     CLEANUP(
  1082.     Vertex, alloc_Vertex, (!ptr->in),
  1083.     volume->firstVertex, volume->lastVertex, volume->numVertex
  1084.     );
  1085.     CLEANUP(
  1086.     Segment, alloc_Segment, (ptr->state==CULLED_OUT),
  1087.     volume->firstSegment, volume->lastSegment, volume->numSegment
  1088.     );
  1089.     CLEANUP(
  1090.     Face, alloc_Face, (ptr->numEdge==0),
  1091.     volume->firstFace, volume->lastFace, volume->numFace
  1092.     );
  1093. }
  1094.  
  1095. /*---------------------------------------------------------------------------*/
  1096. /*    deleteVolume
  1097. /*---------------------------------------------------------------------------*/
  1098. static void
  1099. deleteVolume( Volume *volume )
  1100. {
  1101.     int j,numf;
  1102.     Face *face;
  1103.  
  1104.     if (volume == NULL) return;
  1105.  
  1106. #define    DELETE( Type, Func, numItem, firstItem )\
  1107.     {\
  1108.     int i, num = numItem;\
  1109.     Type *ptr = firstItem;\
  1110.     for (i=0; i<num; i++)\
  1111.     {\
  1112.         Type *next = ptr->next;\
  1113.         Func( ptr );\
  1114.         ptr = next;\
  1115.     }\
  1116.     }
  1117.  
  1118.     DELETE( Vertex, alloc_Vertex, volume->numVertex, volume->firstVertex );
  1119.     DELETE( Segment, alloc_Segment, volume->numSegment, volume->firstSegment );
  1120.  
  1121.     numf = volume->numFace;
  1122.     face = volume->firstFace;
  1123.     for (j=0; j<numf; j++)
  1124.     {
  1125.     Face *nextface = face->next;
  1126.     DELETE( Edge, alloc_Edge, face->numEdge, face->firstEdge );
  1127.     alloc_Face( face );
  1128.     face = nextface;
  1129.     }
  1130.  
  1131.     alloc_Volume( volume );
  1132. }
  1133.  
  1134. /*---------------------------------------------------------------------------*/
  1135. /*    newSegment
  1136. /*---------------------------------------------------------------------------*/
  1137. static Segment *
  1138. newSegment( Volume *volume, Vertex *v0, Vertex *v1 )
  1139. {
  1140.     int i, num = volume->numSegment;
  1141.     Segment *segment = volume->firstSegment;
  1142.  
  1143.     for (i=0; i<num; i++, segment=segment->next)
  1144.     {
  1145.     if (((segment->start==v0) && (segment->end==v1)) ||
  1146.         ((segment->start==v1) && (segment->end==v0)))
  1147.         return segment;
  1148.     }
  1149.  
  1150.     segment = alloc_Segment( NULL );
  1151.     if (!volume->numNewSegment)
  1152.     volume->firstNewSegment = segment;
  1153.     volume->numNewSegment++;
  1154.  
  1155.     segment->state = CULLED_IN;
  1156.     segment->start = v0;
  1157.     segment->end = v1;
  1158.     if (volume->lastSegment == NULL)
  1159.     {
  1160.     volume->firstSegment = volume->lastSegment = segment;
  1161.     segment->last = segment->next = NULL;
  1162.     }
  1163.     else
  1164.     {
  1165.     volume->lastSegment->next = segment;
  1166.     segment->last = volume->lastSegment;
  1167.     volume->lastSegment = segment;
  1168.     }
  1169.     segment->next = NULL;
  1170.     volume->numSegment++;
  1171. }
  1172.  
  1173. /*---------------------------------------------------------------------------*/
  1174. /*    newFace
  1175. /*---------------------------------------------------------------------------*/
  1176. static Face *
  1177. newFace( int numEdge, GLdouble *equation )
  1178. {
  1179.     int i;
  1180.     Face *face = alloc_Face( NULL );
  1181.     Edge *edge;
  1182.  
  1183.     face->numEdge = numEdge;
  1184.     edge = face->firstEdge = alloc_Edge( NULL );
  1185.     for (i=1; i<numEdge; i++)
  1186.     {
  1187.     edge->next = alloc_Edge( NULL );
  1188.     edge->next->last = edge;
  1189.     edge = edge->next;
  1190.     }
  1191.     edge->next = face->firstEdge;    /* ring buffer */
  1192.     face->firstEdge->last = edge;
  1193.  
  1194.     if (equation != NULL)
  1195.     for (i=0; i<4; i++)
  1196.         face->equation[i] = equation[i];
  1197.  
  1198.     return face;
  1199. }
  1200.  
  1201. /*---------------------------------------------------------------------------*/
  1202. /*    getMinMaxLat
  1203. /*---------------------------------------------------------------------------*/
  1204. /*
  1205.     get min and max z-values of volume
  1206. */
  1207. static int
  1208. getMinMaxLat( Volume *volume, float *latmin, float *latmax )
  1209. {
  1210.     int i, num;
  1211.     Vertex *v;
  1212.     float min, max;
  1213.  
  1214.     if (volume == NULL) return 0;
  1215.     if ((num=volume->numVertex) <= 0) return 0;
  1216.     if ((v=volume->firstVertex) == NULL) return 0;
  1217.  
  1218.     min = max = v->z;
  1219.     v = v->next;
  1220.     for (i=1; i<num; i++, v=v->next)
  1221.     {
  1222.     min = INF(min, v->z );
  1223.     max = SUP(max, v->z );
  1224.     }
  1225.  
  1226.     *latmin = min;
  1227.     *latmax = max;
  1228. }
  1229.  
  1230. /*---------------------------------------------------------------------------*/
  1231. /*    getMinMaxLon
  1232. /*---------------------------------------------------------------------------*/
  1233. /*
  1234.     for a given z=cste line, get min&max x-values of volume
  1235. */
  1236. static int
  1237. getMinMaxLon( Volume *volume, float lat, float *lonmin, float *lonmax )
  1238. {
  1239.     int    i, num, first_found = 0;
  1240.     float    dl, lon, lonA, latA, lonB, latB, max, min;
  1241.     Segment    *segment;
  1242.  
  1243.     if (volume == NULL) return 0;
  1244.     if ((num=volume->numSegment) <= 0) return 0;
  1245.     if ((segment=volume->firstSegment) == NULL) return 0;
  1246.  
  1247.     for (i=0; i<num; i++, segment=segment->next)
  1248.     {
  1249.     lonA = segment->start->x;
  1250.     latA = segment->start->z;
  1251.     lonB = segment->end->x;
  1252.     latB = segment->end->z;
  1253.     
  1254.     if ((lat==latA) && (latA==latB))
  1255.     {
  1256.         if (first_found)
  1257.         {
  1258.         min = INF( min, INF(lonA,lonB));
  1259.         max = SUP( max, SUP(lonA,lonB));
  1260.         }
  1261.         else
  1262.         {
  1263.         min = INF(lonA,lonB);
  1264.         max = SUP(lonA,lonB);
  1265.         first_found = 1;
  1266.         }
  1267.     }
  1268.     else if (latA != latB)
  1269.     {
  1270.         dl = latB - latA;
  1271.         if (((dl>0) && (lat>=latA) && (lat<=latB)) ||
  1272.         ((dl<0) && (lat<=latA) && (lat>=latB)))
  1273.         {
  1274.         lon = lonA + (lonB-lonA)*((lat-latA)/dl);
  1275.         if (first_found)
  1276.         {
  1277.             max = SUP( max, lon );
  1278.             min = INF( min, lon );
  1279.         }
  1280.         else
  1281.         {
  1282.             max = min = lon;
  1283.             first_found = 1;
  1284.         }
  1285.         }
  1286.     }
  1287.     }
  1288.  
  1289.     *lonmin = min;
  1290.     *lonmax = max;
  1291.     return first_found;
  1292. }
  1293.  
  1294. #define CHUNKSIZE_VOL 1
  1295.  
  1296. Volume*
  1297. alloc_Volume(
  1298.     Volume *ptr
  1299.     ) {
  1300.     static Volume *store;
  1301.     static int free=0;
  1302.  
  1303.     if(!ptr) {
  1304.     if(--free < 0) {
  1305.         register int i;
  1306.         if((store = (Volume*)malloc(CHUNKSIZE_VOL*sizeof(Volume))) == NULL)
  1307.         {fprintf(stderr, "SYSTEM ERROR: malloc (exiting)"); exit(1);}
  1308.         for(i=0, ptr=store; i<CHUNKSIZE_VOL-1; ++i, ++ptr)
  1309.         ptr->next = ptr+1;
  1310.         free += CHUNKSIZE_VOL;
  1311.         }
  1312.     ptr = store;
  1313.     store = store->next;
  1314.     }
  1315.     else {ptr->next = store; store = ptr; free++;}
  1316.     return(ptr);
  1317. }
  1318.  
  1319. #define CHUNKSIZE_VERT 12
  1320.  
  1321. Vertex*
  1322. alloc_Vertex(
  1323.     Vertex *ptr
  1324.     ) {
  1325.     static Vertex *store;
  1326.     static int free=0;
  1327.  
  1328.     if(!ptr) {
  1329.     if(--free < 0) {
  1330.         register int i;
  1331.         if((store = (Vertex*)malloc(CHUNKSIZE_VERT*sizeof(Vertex))) == NULL)
  1332.         {fprintf(stderr, "SYSTEM ERROR: malloc (exiting)"); exit(1);}
  1333.         for(i=0, ptr=store; i<CHUNKSIZE_VERT-1; ++i, ++ptr)
  1334.         ptr->next = ptr+1;
  1335.         free += CHUNKSIZE_VERT;
  1336.         }
  1337.     ptr = store;
  1338.     store = store->next;
  1339.     }
  1340.     else {ptr->next = store; store = ptr; free++;}
  1341.     return(ptr);
  1342. }
  1343.  
  1344. #define CHUNKSIZE_EDGE 25
  1345.  
  1346. Edge*
  1347. alloc_Edge(
  1348.     Edge *ptr
  1349.     ) {
  1350.     static Edge *store;
  1351.     static int free=0;
  1352.  
  1353.     if(!ptr) {
  1354.     if(--free < 0) {
  1355.         register int i;
  1356.         if((store = (Edge*)malloc(CHUNKSIZE_EDGE*sizeof(Edge))) == NULL)
  1357.         {fprintf(stderr, "SYSTEM ERROR: malloc (exiting)"); exit(1);}
  1358.         for(i=0, ptr=store; i<CHUNKSIZE_EDGE-1; ++i, ++ptr)
  1359.         ptr->next = ptr+1;
  1360.         free += CHUNKSIZE_EDGE;
  1361.         }
  1362.     ptr = store;
  1363.     store = store->next;
  1364.     }
  1365.     else {ptr->next = store; store = ptr; free++;}
  1366.     return(ptr);
  1367. }
  1368.  
  1369. #define CHUNKSIZE_FACE 8
  1370.  
  1371. Face*
  1372. alloc_Face(
  1373.     Face *ptr
  1374.     ) {
  1375.     static Face *store;
  1376.     static int free=0;
  1377.  
  1378.     if(!ptr) {
  1379.     if(--free < 0) {
  1380.         register int i;
  1381.         if((store = (Face*)malloc(CHUNKSIZE_FACE*sizeof(Face))) == NULL)
  1382.         {fprintf(stderr, "SYSTEM ERROR: malloc (exiting)"); exit(1);}
  1383.         for(i=0, ptr=store; i<CHUNKSIZE_FACE-1; ++i, ++ptr)
  1384.         ptr->next = ptr+1;
  1385.         free += CHUNKSIZE_FACE;
  1386.         }
  1387.     ptr = store;
  1388.     store = store->next;
  1389.     }
  1390.     else {ptr->next = store; store = ptr; free++;}
  1391.     return(ptr);
  1392. }
  1393.  
  1394. #define CHUNKSIZE_SEGMENT 16
  1395.  
  1396. Segment*
  1397. alloc_Segment(
  1398.     Segment *ptr
  1399.     ) {
  1400.     static Segment *store;
  1401.     static int free=0;
  1402.  
  1403.     if(!ptr) {
  1404.     if(--free < 0) {
  1405.         register int i;
  1406.         if((store = (Segment*)malloc(CHUNKSIZE_SEGMENT*sizeof(Segment))) == NULL)
  1407.         {fprintf(stderr, "SYSTEM ERROR: malloc (exiting)"); exit(1);}
  1408.         for(i=0, ptr=store; i<CHUNKSIZE_SEGMENT-1; ++i, ++ptr)
  1409.         ptr->next = ptr+1;
  1410.         free += CHUNKSIZE_SEGMENT;
  1411.         }
  1412.     ptr = store;
  1413.     store = store->next;
  1414.     }
  1415.     else {ptr->next = store; store = ptr; free++;}
  1416.     return(ptr);
  1417. }
  1418.